arXiv:1502.02869vl [physics.flu-dyn] 10 Feb 2015 


Under consideration for publication in J. Fluid Mech. 


1 


Universal mechanism for air entrainment 
during liquid impact 

Maurice H.W. Hendrix^’^, Wilco Bouwhuis^, Devaraj van der Meer^, 
Detlef Lohse^, and Jacco H. Snoeijer^’^ 

^ Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AE 

Enschede, The Netherlands, 

^ Laboratory for Aero and Hydrodynamics, Delft University of Technology, Leeghwaterstraat 
21, NL-2628 CA Delft, The Netherlands, 

® Mesoscopic Transport Phenomena, Eindhoven University of Technology, Den Dolech 2, 5612 

AZ Eindhoven, The Netherlands 

(Received 11 February 2015) 

When a mm-sized liquid drop approaches a deep liquid pool, both the interface of the 
drop and the pool deform before the drop touches the pool. The build up of air pressure 
prior to coalescence is responsible for this deformation. Due to this deformation, air can 
be entrained at the bottom of the drop during the impact. We quantify the amount of 
entrained air numerically, using the Boundary Integral Method (BIM) for potential flow 
for the drop and the pool, coupled to viscous lubrication theory for the air film that has 
to be squeezed out during impact. We compare our results to various experimental data 
and find excellent agreement for the amount of air that is entrapped during impact onto 
a pool. Next, the impact of a rigid sphere onto a pool is numerically investigated and the 
air that is entrapped in this case also matches with available experimental data. In both 
cases of drop and sphere impact onto a pool the numerical air bubble volume 14 is found 
to be in agreement with the theoretical scaling VblVdrop/sphere ^ where St is the 

Stokes number. This is the same scaling that has been found for drop impact onto a solid 
surface in previous research. This implies a universal mechanism for air entrainment for 
these different impact scenarios, which has been suggested in recent experimental work, 
but is now further elucidated with numerical results. 


1. Introduction 

The impact of a drop or a solid sphere onto a liquid pool can encompass various types 
of air entrainment. One possibility is that air is entrained at the top of the impacting 
object when the crater that is created during impact collapses, see for example Oguz & 
Prosperetti (1990); Pumphrey & Elmore (1990); Wang et al. (2013); Chen & Guo (2014). 
Another type of air entrainment may occur at the bottom of the impacting object: the 
thin air film that is squeezed out at the impact zone is accompanied by a pressure in¬ 
crease that deforms the interface of the liquid before the impacting object touches the 
pool, which may result in air entrapment (Marston et al. 2011; Tran et al. 2013). The 
early stages of deformations can be described analytically (Bouwhuis et al. -). In case the 
impacting object is a drop, instead of a single entrapped bubble, also a collection of micro¬ 
scopic bubbles may be entrapped which can create intriguing morphologies (Thoroddsen 
et al. 2012). This is also referred to as Mesler entrainment (Esmailizadeh & Mesler 1986; 
Pumphrey & Elmore 1990). The same mechanism that is responsible for bubble entrap¬ 
ment at the bottom of an impacting object on a pool holds for of air entrapment at the 
bottom of an impacting drop onto a solid (van Dam & Le Clerc 2004; Mani et al. 2010; 
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Figure 1: Air bubble entrapment for different impact scenarios. Bubbles and deformations 
are not drawn to scale, (o) Rigid sphere impact onto a pool. The pool deforms due to an 
increase in air pressure right under the sphere before it touches the pool, which results 
in an entrapped air bubble, {h) Drop impact onto a pool. Not only the pool, but also the 
drop consists of a deformable interface. As a result, the increased air pressure deforms 
both the pool and the drop and an air bubble is entrapped, (c) Drop impact onto a solid. 
Also here, a local increase in air pressure deforms the drop before it touches the solid 
and results in an entrapped air bubble. 


Hicks & Purvis 2010; Bouwhuis et al. 2012). In fact, the initial geometry of the problems 
is identical, see figure 1 in which the different impact scenarios and air entrapment have 
been depicted. We also refer to figure 5 of Tran et al. (2013), who first worked out this 
analogy. 

Previously, air bubble entrapment for drop impact onto a solid surface has been quanti¬ 
fied experimentally, theoretically, and numerically (Mandre et al. 2009; Mani et al. 2010; 
Hicks & Purvis 2010, 2011; Bouwhuis et al. 2012). If the effect of surface tension can be 
neglected, the following scaling for the entrapped air bubble volume was found: 


Vb/Vdrop - St-^/\ ( 1 . 1 ) 

Here VbjVdrop is the air bubble volume normalized by the drop volume and St is the 
Stokes number which is defined as St = piRU/rjg, where pi is the liquid density, R the 
droplet radius, U its impact velocity, and ijg is the viscosity of the surrounding gas, in 
this case air. The Stokes number represents the competing effect of the viscous force of 
the draining air film and the inertial force of the liquid which ultimately determines the 
air bubble volume. The same scaling was found experimentally for impact of a sphere 
onto a pool (Marston et al. 2011), and a drop onto a pool (Tran et al. 2013). 

In this paper we try to capture the mechanism of air entrapment during impact onto 
a pool numerically. We will employ a boundary integral method (BIM) for potential 
flow describing the liquid phase coupled to viscous lubrication theory for the draining 
microscopic air film. The advantage of using a boundary integral method becomes evident 
when the interface of the impacting object comes close to the pool and one has to resolve 
the microscopic air layer together with the macroscopic liquid scale. This difference in 
length scale can be a thousandfold for the case of a millimeter sized drop impacting 
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onto a pool squeezing out an air film with a typical thickness of micrometer. In fact, the 
difference in length scale in the final stages of impact diverges to infinity as the drop is 
about to coalesce with the pool. Using a boundary integral method guarantees excellent 
interface representation, since all variables such as liquid velocity and pressure are defined 
at the interface. At the same time, the computational cost is modest, since the boundary 
integral method allows the potential problem to be solved only at the boundaries of the 
liquid domain: quantities in interior points can be calculated optionally as a function 
of the solution at the boundary. To achieve the same accurate interface representation 
and solving the full Navier-Stokes equations, using for example a volume-of-fluid method 
(see for example Thoraval et al. (2012); Guo et al. (2014)), would require a much larger 
computational cost. 

In section 2 we explain the theoretical framework together with the numerical method. 
In section 3 we will present the results of numerical simulation: we will identify details 
of the pressure development in the air film and deformation of the interfaces at the 
impact zone. The results of the numerical model will be compared with available results 
regarding the entrapped bubble volume from multiple experimental works and will be 
compared with the scaling law equation (1.1). We conclude with section 4 in which also 
suggestions for further research are discussed. 


2. Theory 

2.1. Dimensional analysis and numerical method 

The Reynolds number of the liquid drops we model, which is defined as Rei = piRU/rji, 
is assumed to be large, Rei ^ 1. Here pi and ry are respectively the density and the 
dynamic viscosity of the liquid, U is the impact velocity and R is the radius of drop. The 
flow can thus be regarded irrotational, that is V x u = 0. Under the additional constraint 
of incompressible flow inside the drop this allows the liquid dynamics to be modeled with 
a harmonic function </>, to which the velocity field u is related through: 


u = V(j) ( 2 . 1 ) 

The fact that the velocity potential (f obeys the Laplace equation V^(/) = 0 is used to 
efficiently solve the potential problem, and thus the dynamics of the liquid, using the 
Boundary Integral Method (BIM). While the Reynolds number of the drop is large, 
the Reynolds number of the thin gaseous air layer Reg = PgHdU/ijg is typically small, 
Reg <C 1. Here pg is gas density and Hd is the air film thickness in the center of the film 
which is referred to as the dimple height. The length scale characterizing the air layer in 
the lateral extension of the air film is denoted by L, see figure 2a. As shown in Bouwhuis 
et al. (2012), Hd ^ L which in combination with the low Reynolds number of the gas 
allows the film to be described with viscous lubrication theory, see for example Leal 
(1992). The dimensionless group reflecting the presence of air is the Stokes number St = 
piRU/rjg which compares the viscous force of the air layer to the inertial force in the 
drop. This number is relevant for describing dimple formation, since, for high enough 
impact velocity U, this process is determined by two competing forces: the force of the 
viscous air layer trying to deform the drop in the center and the opposing inertial force 
of the drop, which must be slowed down locally in order to form a dimple. Additional 
dimensionless numbers incorporating surface tension 7 are the Weber number We and 
the capillary number Ca based on the gas properties. Summarizing, we thus have the 
following dimensionless parameters: 
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c.l? 

m Vg Vg 7 St 

( 2 . 2 ) 

The impact of a liquid drop onto a pool of the same liquid and the impact of a rigid 
sphere onto a liquid pool can be described with the same dimensionless numbers. As the 
initial geometry of the problems is identical, the difference lies in the deformability of 
the object, which is zero in case of the solid. The two effective control parameters that 
we will use here in our theoretical framework are St and We. 

In figure 2a an illustration of the impact of a drop onto a pool, together with the used 
method is shown. As is clear from this figure, the coupling between the dynamics of the 
air layer and the dynamics of the liquid is essential since the two liquid domains feel each 
other through the pressure build up in the viscous air layer. The lubrication pressure 
Pg acts on the liquid surface and appears in the unsteady Bernoulli equation which 
serves as a boundary condition in the BIM which is applied at the liquid surface, see 
also Bouwhuis et al. (2012). As we have two liquid domains, two separate BI equations 
are solved. We take the width of the pool large enough to approach the dynamics of 
an infinite liquid pool. In this case a width of 4.5 times the drop radius was found to 
be sufficient. We focus on quantifying the amount of entrapped air by integrating the 
enclosed air pocket up to the moment the air layer reaches a physical minimum thickness 
of 0.4 |xm. At this point the volume of the enclosed air has converged and a subsequent 
rupture of the air film will prevent further drainage which results in an entrapped air 
bubble (Bouwhuis et al. 2012). As we focus on the dynamics just prior to rupture we can 
make use of an axisymmetric framework. We restrict ourselves to the inertial regime for 
which experimental results (Tran et al. 2013) are available for a direct comparison. Since 
the air layer continually deforms and translates during the impact, lubrication equations 
have been developed in a moving coordinate system which is aligned with the interface 
of the drop. These equations will be derived in the next section. 

2.2. Lubrication in moving and tilted coordinate system 

In this section we develop an expression for the pressure Pg in the air film based on 
lubrication theory in a moving (||,T)-coordinate system which is aligned along the drop 
surface, see the sketch in figure 2. The reason for doing this (rather than just using 
the standard (r,z)-coordinate system) is that especially for the drop onto pool impact, 
the moving (||,T)-coordinate system is not necessarily oriented as the (r, 2 )-coordinate 
system and therefore only the first guarantees an accurate description of the draining air 
film. The drop surface is taken as a reference, and the curvilinear coordinate || is defined 
along the drop, starting at the axis of symmetry. At some large radial coordinate ||oo 
we assume atmospheric pressure. The coordinate perpendicular to || is defined to be T. 
The gap height h(r, t) is defined as the length of the perpendicular line from the drop 
projected onto the liquid pool. The two surfaces in the impact zone are assumed to be 
nearly parallel (|9||/i| <C 1), so we can apply lubrication theory. 

It can be shown (see appendix A) that the continuity equation in this new (||,T)- 
coordinate system reads: 


— + i9||U|l + = 0. (2.3) 

At the interface of the liquid pool (T= h) we know that the fluid particles have to 
move with the interface. This is mathematically described with the kinematic boundary 
condition: 
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Figure 2: (a) Schematic of drop impact onto a a pool. The used methods are indicated 
in the figure: both the liquid domains are modelled with potential flow, while the air 
layer is described with Stokes flow. The gray arrows indicate that the flow of the air 
film is coupled to the dynamics of the liquid domains and vice versa. (6) Definition of 
the (||,_L)-coordinate system where || is aligned along the drop curve and T is the unit 
normal with respect to the drop. 


dth+ (m|| -u_L|x=o- (2-4) 

Here dth is the time derivative of h. We now integrate equation (2.3) along the gap height 
h and obtain: 


nh ^ ph ph 

J -yd±+J (9||U||dT=-y aj_u_Ld T= m_l|_l=o 

Using Leibniz integral rule for the second integral on the left hand side we find: 


(2.5) 


Ur 

/ — dT+9||/ M|| d T - (m|| a||/i) = ■u_L|_L=o - u±|±=;i. (2.6) 

Jo ^ Jo 

We now use the kinematic boundary condition formulated in equation (2.4) for the third 
term on the LHS to obtain: 


Ur 

/ —dT+i9||/ ui\d±+dth-u±i±=h + u±i±=o = u±i±=o-u±i 

Jo ’’’ Jo 

Canceling the terms and on both sides gives: 


_L=/f 


(2.7) 


nh ph 

yd T+5,1 


M|l d T d-dth = 0. (2.8) 

/o ' Jo 

We still have to describe Ur within the new (||,_L)-coordinate system. Therefore we sub¬ 
stitute Ur = Mj_ cos 9 — u\\ sin 9 in the equation above to get: 


1 


1 


-M_LCOS0dT— / -Mil sin^d T-|-9| 


I / Mil d T -\-dth = 0. (2.9) 

Jo ' ^0 ' Jo 

We assume that the main flow of the air that is squeezed out from the gap is along 
the II coordinate, which implies that m_l is relatively small, so we neglect the first term. 
The second term is an integral with respect to T containing the variable r. This radial 
coordinate r across /i is a function of T: r =T cos0 -I- c(||). Here c(||) is the value of r at 
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the drop surface (_L= 0) for some coordinate ||. We thus substitute this expression for r 
into equation (2.9) and neglect the first term to find: 

d _L +dth = 0. (2.10) 

2.2.1. Flow profile within the air film 

As has been previously described the Reynolds number of the thin air film is small, 
Rcg <C 1, and the geometry of the problem, Hd L, allows us to use lubrication theory. 
In the (||,_L)-coordinate system, the Stokes equations can then be reduced to: 

d\\Pg =Vgdiu\\. (2.11) 

We can integrate equation equation (2.11) twice with respect to _L, employing a no slip 
boundary condition at the drop surface (m|| = Ud at _L= 0) as well as at the surface of 
the pool (mii = Up at _L= h): 

«|| = (^([/p - Ud)j + Ud^ + ^ h). (2.12) 

The first term of equation (2.12) can be associated with Couette flow, caused by the 
movement of the interfaces. The second term can be associated with Poiseuille flow, which 
is driven by the radial pressure gradient, see also Klaseboer et al. (2000). Substituting 
this expression for My in our equation for mass conservation, equation (2.10), we get: 


sind 


/q T cos 9 + c{ 


-Mil d T -|-9|| 


sin 9 


T cos 9 + c 


{Up - Ud)— + Ud] + - —9||(T^ - Fh) 

■^Vg 


d T 



T 


1 


Ud)T + Ud]+^^d^\Pg{P 


T h) 


d T +dth = 0. 


(2.13) 


In the first integral we deal with a prefactor sin 9/(1. cos 9+c). When taking the geometry 
of the problem into account we note that T cosd <C c. We can thus write sin0/(T 
cos 9 -I- c) w sm9/c. Performing the integrals of equation (2.13) under this assumption 
yields: 


_s^ (Up + Ud) - + ^11 (Up + Ud) - + 9th - 0. (2.14) 

If we define Gdl) = {Up + Ud) — we can transform the above equation into 

a first order inhomogeneous linear ODE for G(||): 


G(||)-a(||)G(||) = /(||). 
Here a(||) and /(||) are known functions of ||: 


«(ll) = 


sind 


(2.15) 


(2.16) 


/(II) = -dth 


( 2 . 17 ) 
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2.2.2. Solving the first order inhomogeneous ODE for G(||) 

Equation (2.15) can be solved with help of an integrating factor I defined as /(||) = 
e-/“(ll)<i|l. Using the boundary condition G(||) = 0 for ||= 0, because we have zero 
pressure gradient in the center of symmetry, and also zero tangential velocities, we can 
multiply equation (2.15) with /(||) and solve for G(||): 


G(ii) = y^(^"/(y)/(y)dy) (2.18) 

with /(II) = e“ -^0 '^11. 

We can now substitute G(||) = {Up + Ud) — T|^'9||/’g^ back into equation (2.18) to 
obtain an equation for d\\Pg\ 


d\\P9 = ^(i)/(y) di j + (2-19) 

We note that we have to evaluate two numerical integrals to calculate di\Pg. In order 
to find the pressure Pg{\\) we integrate equation (2.19) using atmospheric pressure for 
some large value for ||oo well outside the thin air gap as a boundary value. As a check 
of our analysis we now orientate the (||,_L)-coordinate system in such away that ||= r, 
to recover the lubrication equation in the conventional (r, 2 ;)-coordinate system. In that 
case we have 0 = —7r/2, and we can write for a(||): 


a(||= r) 


The integrating factor I now becomes: 


sin0 

r 


I 

r 


( 2 . 20 ) 


/(II = r) = e-dJ'AUd? ^ ginr ^ ^ (2.21) 

Substituting equation (2.21) into equation (2.19) and using the proposition ||= r and 
setting Ub = 0 and Ud = 0, we can now write equation (2.19) as: 


drPg = - 


h^ 


iw) (L - 5 <^* + = ^ (; 


rdfhdr 


( 2 . 22 ) 

We inspect that this equation (2.22) is the equation for the radial pressure gradient for 
viscous lubrication theory in the conventional (r, 2 )-coordinate system (Bouwhuis et al. 
2013), which gives a consistency check for our analysis. This was also numerically verihed. 


3. Results 

In this section simulation results will be discussed, starting with section 3.1 in which 
the drop impact onto a pool will be treated. The interface deformations and pressure 
development in the viscous air layer will be quantihed. In section 3.2 we will focus on 
rigid sphere impact onto a pool. For both impact scenarios we will quantify the size of the 
air bubble that is entrapped and directly compare with various experimental results (Tran 
et al. 2013; Marston et al. 2011; Bouwhuis et al. 2012). In section 3.3 we will compare 
the dynamics of both impact scenarios and identify symmetrical behavior. 
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3.1. Drop impact onto a pool 

Figure 3a displays a typical result for drop impact onto a pool. The results are expressed 
in dimensional form to match the experimental conditions of the work of Tran et al. 
(2013), to which the numerical results in this work will be compared. In the hrst frame 
corresponding to t = 0 ms the initial condition of the simulation at the impact zone 
is shown. An initial separation of Iiq = 50 pm is used. Convergence tests regarding the 
initial release height have been conducted, and an initial separation of ho = 50 pm was 
found to be appropriate for the lubrication pressure to be still negligible at this distance 
for the parameter range which is of interest in this study. At t = 0.12 ms it can be seen 
that the pool and the drop experience the increased air pressure and thus the interfaces 
deform. In the lower panel of this frame, the increase in pressure is indeed visible. At 
t = 0.15 ms the drop is getting closer to the pool, and the interfaces have been further 
deformed. It can also be noted that the pressure maximum corresponds to a location 
where the separation between the drop and the pool is smallest. The location of smallest 
separation is now not located in the center at r = 0 anymore. This behavior is typical 
for impact events involving a free surface and has been experimentally observed for e.g. 
drop impact onto a pool (Thoroddsen et al. 2012; Tran et al. 2013), drop impact onto 
a solid surface (Van der Veen et al. 2012; De Ruiter et al. 2012), sphere impact onto a 
pool (Marston et al. 2011) and bubble impact onto a wall in a liquid tank (Hendrix et al. 
2012). In the final frame t = 0.17 ms we observe that the two interfaces are very close 
together having a minimum separation of 0.4 pm. 

We note that the interfaces up to the final stage of impact are very well resolved, 
see hgure 3b in which the final frame at t = 0.17 ms is shown on various scales while 
keeping both axes the same length scale. In the first frame in this figure a macroscopic 
view of the simulation domain is shown. In the second frame the impact zone is selected 
and magnified. The slender geometry of the microscopic air film can be noted. In the 
third frame the region of closest separation is magnified. Indeed, the interfaces are very 
close together, the minimum separation is 0.4 pm. The computational nodes used for 
discretization of the surface are also shown in this final frame. An adaptive grid on the 
fluid surface allows for local refinement at the region of closest separation which results in 
the total number of nodes to be only of order 100, while capturing both the microscopic 
dynamics at the impact zone and the large scale motion of the millimeter sized drop. 
Note that the slanted orientation of the free surfaces in figure 3b justifies the need of 
using a description in terms of the (||,T)- rather than the (r,z)-coordinate system. 

We further note from the hnal frame in figure 3a that a microscopic air film hnds itself 
trapped between the drop and the pool. It is this entrapped air that constitutes the air 
bubble that is dragged into the liquid when the air film ruptures at the thinnest point and 
breaks the axisymmetry of the problem. In this work we do not attempt to simulate the 
complex rupture process of the air film itself, which is ultimately determined by surface 
chemistry, see for example Saylor & Bounds (2012) who point out the importance of 
cleanliness of the fluid surface with respect to film rupture. Instead, we focus on the 
dynamics up to the rupture point, which is taken to happen at a rupture height of 
0.4 pm. At this point, the volume of the entrained air has converged and can thus be 
determined, see the final frame of figure 3a. This procedure is in line with previous 
research (Bouwhuis et al. 2012), where experimentally the volume of the air pocket just 
before rupture was indeed found identical to the volume of the entrapped bubble. 

3.2. Rigid sphere impact onto a pool 

The impact of a sphere onto a pool prior to coalescence is similar to the case of a drop 
impacting on a pool, except that in the case of an impacting sphere the deformability of 
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t = 0.00 ms t = 0.13 ms t = 0.15 ms 1 = 0.17 ms 



r (mm) r (mm) r (mm) r (mm) 






Figure 3: (a) Drop impact onto a liquid pool. Note the different length scales for the 
r-axis and z-axis in the shape plots. The impact speed is 17 = 0.42 m/s and the drop 
radius is i? = 0.95 mm. The density and surface tension of the liquid are respectively p = 
916 kg/m^ and 7 = 0.020 N/m. These impact parameters correspond to St = 2.0 x 10^ 
and We = 7.7. The simulation starts at time t = 0 ms at a separation of hr=o = 50 pm. 
Due to the approach of the sphere, the excess air pressure Pg will increase and acts 
on both the drop and the liquid pool (t = 0.13 ms). At the final stage {t = 0.17 ms) 
the minimum separation of the interfaces reaches 0.4 pm and the simulation is stopped. 
The bubble volume Vf, can thus be determined, (b) Part of the simulation domain with 
detailed snapshots of the air film at t = 0.17 ms. In the third snapshot the actual node 
distribution around the smallest separation point can be inspected. This is the most 
refined distribution of computational nodes that is used. For the region outside the gap 
a coarser node distribution is sufficient. 
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Figure 4: Rigid sphere impact onto a pool. The impact speed is U = 0.42 m/s and the 
radius is R = 0.95 mm. The density and surface tension of the fluid are respectively p = 
916 kg/m^ and 7 = 0.020 N/m. These impact parameters correspond to St = 2.0 x 10^ 
and We = 7.7. 


the impacting object is zero. This scenario has been simulated by letting an undeformable 
sphere approach the pool. The same equations are solved as described in Section 2, except 
no BIM is needed for the impacting sphere since the interface of the sphere is fixed. The 
result is depicted in figure 4. Just as in the case of drop impact onto a pool, a microscopic 
air bubble is entrapped. As can be inspected, the air bubble has a similar shape, but 
its size is smaller than in case a drop impacts onto the pool, as can be inferred from a 
comparison to figure 3a. 

The size of the air bubble can be quantified from the numerical simulation and is 
compared for both drop impact and sphere impact onto a pool with various experimental 
results in figure 5a. We see that the numerical results of both drop and sphere impact 
onto a pool are in quantitative agreement with experimental work. We note that the 
numerical results show that the air bubble volume is indeed larger when a drop instead 
of a sphere impacts onto a pool for all St, which is supported by experiments of Tran 
et al. (2013). Furthermore, we observe that numerical results are in agreement with the 
scaling law presented in equation (1.1), VbjVdrop St~'^^^. As experiments have shown, 
in this regime, viscosity of the liquid is not important for the final bubble volume that is 
entrapped, see Marston et al. (2011); Tran et al. (2013), which is again confirmed by the 
current modeling technique which captures the essential physics which determine the air 
bubble volume: a potential flow calculation that does not involve liquid viscosity coupled 
to viscous lubrication theory for the intervening airlayer. 

3.3. Deformations of interfaces: symmetrical behavior 
We will now further investigate the fact that the bubble volume for drop impact onto 
a pool is larger compared to the case where we deal with only one deformable interface 
during impact as is the case with rigid sphere impact onto a pool. In figure 6 a closer 
inspection of the drop impact onto a pool is depicted. In this figure we track the relative 
deformation of both the pool and the drop, denoted by Sdrop and Spool respectively. 
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BIM, rigid sphere on pooi 

15mm steel sphere on water (Marston et al., 2011) 
20mm steel sphere on 10cSt oil (Marston et al., 2011) 
20mm steel sphere on 25cSt oil (Marston et al., 2011) 
20mm steel sphere on IcSt SDS (Marston et al., 2011) 
Ethanol drop on solid surface (Bouwhuis et al., 2012) 


• BIM, drop on pool 

• Drop on pool, silicone oil (5cSt) (Tran et al. 2013) 

• Drop on pool, silicone oil (lOcSt) (Tran et al. 2013) 
■ Drop on pool, silicone oil (20cSt) (Tran et al. 2013) 


Figure 5: Figure adapted from al. Tran et al. (2013). BIM results are superimposed 
in yellow symbols, (a) Various experimental data for the normalized bubble volume 
ybiVsphere/drop sre shown. Excellent quantitative agreement was found with numerical 
results, (b) The data, both numerical and experimental, was found to collapse on one 
single curve by normalizing 14 as Vb j{nVsphere/drop) with n the number of free interfaces 
involved during impact. This number is 2 instead of 1 in case of drop impact onto a liquid 
pool. 


Here 5drop is defined as the deformation of the drop relative to an undeformed sphere 
impacting with constant speed U and Spool is defined as the deformation of the pool 
relative to horizon z = 0. Interestingly, we note that both interface deformations behave 
identically. One may expect that two deformable interfaces which react similar to an 
external pressure, deform in an identical way. Note, however, that the upper domain 
(drop) and the lower domain (pool) do not have the same unperturbed geometry, owing 
to the radius of curvature of the drop. Since both media respond identically to the 
pressure pulse, the weak curvature with respect to the width of the localized pressure 
has a negligible influence: on the scale of the pressure pulse, both domains are essentially 
flat. We therefore expect to recover a symmetric response in the upper and lower domains. 
To illustrate this further, we compute the kinetic energy and the velocity inside the drop 
and the pool using a technique described in Sun et al. (2014) to evaluate quantities close 
to the interface which need special attention as the singular behavior of the Green’s 
function in the Boundary Integral equation becomes apparent for these points. Figure 7a 
shows the result in the frame of the pool. To highlight the symmetry, we also evaluate 
these quantities in a frame moving at a speed U /2 in an upward direction, which results 
in a frame of reference in which both the drop and pool move with a speed U /2 towards 
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Figure 6 : Drop impact onto a pool with a corresponding plot of the relative deformation 
5 of both the pool and the drop. In the final frame of the upper panel the definition 
of <5 is shown. We observe excellent overlap between the relative deformations, which is 
emphasized in the lower panel where A = 6drop — dpooh the difference between the two 
relative deformations, is shown. The same impact conditions as for the case described in 
figure 3 are used: The impact speed is U = 0.42 m/s and the radius is i? = 0.95 mm. 
The density and surface tension of the liquid are respectively p = 916 kg/m^ and 7 = 
0.020 N/m, which corresponds to St = 2.0 x 10^ and We = 7.7. 


each other. Indeed, the velocity fields and kinetic energies are now identically distributed, 
see figure 7b. 

This implies that there will be a bigger entrapped air bubble as compared to the case 
where only one of the interfaces is able to deform. To quantify this hypothesis we compare 
the bubble sizes of drop and sphere impact onto a pool and find a factor 2 difference, 
see figure 5b. Here half the air bubble volume of drop impact onto a pool was found to 
collapse onto the experimental and numerical results incorporating only one deformable 
interface, i.e. the sphere impact onto a pool but also drop impact onto a solid. Tran et al. 
(2013) took another approach to collapse the data of bubble volumes of drop impact 
onto a pool by correcting the corresponding impact St number by a factor 2, which also 
collapses the data. In this present work it is shown that an approach based on considering 
the number of deformable interfaces (either 1 or 2 ) can also serve to obtain a unifying 
view on the air bubble entrapment. 
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Figure 7: Kinetic energy monitoring during drop impact onto a pool, with impact pa¬ 
rameters as described in figure 3. (a) In the left half, the kinetic energy K is color coded: 
black is zero kinetic energy, white is maximum kinetic energy in the system which is 
K = ^pU^. (b) The kinetic energy is recalculated in a moving reference frame moving 
upwards at ^{7. This results in a frame of reference in which both the pool and drop 
move with a speed of ^17 toward each other. Again the left half of the figure shows the 
kinetic energy. We observe a symmetric behavior which supports the hypothesis that the 
pool and drop react in a symmetric way to the local pressure increase. 


4. Conclusion 

In this work air entrapment during liquid drop and rigid sphere impact onto a deep 
liquid pool has been numerically investigated using a Boundary Integral Method (BIM) 
for potential flow for the liquid phase coupled to the viscous lubrication approximation 
for the subphase air which is squeezed out during impact. Excellent agreement with 
experimental work was found when comparing the amount of air that is entrained during 
impact. When considering drop impact onto a pool both liquid interfaces were found 
to deform identically relative to their undeformed shape. This leads to an explanation 
as to why bubble volumes in case of drop impact onto a pool were found to be twice 
the size of those that are found, both experimentally and numerically, in impacts events 
involving only one deformable interface, that is, rigid sphere impact onto a pool and drop 
impact onto a solid. In this study compressibility effects of the air have been neglected. 
It can be expected that at higher impact velocity compressibility of the intervening 
air will be important, see for example Hicks & Purvis (2011). In addition, the current 
modelling technique is limited to an axisymmetric 2D framework. To account for 3D 
impact problems, for which experimental data starts to emerge (Van der Veen et al. 
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2014), the modelling technique needs to be extended to 3D. With a 3D model also 
oblique collisions can be investigated. 


Appendix A. Continuity in curvilinear coordinates 

To derive equation 2.3 in a (||,_L)-coordinate system that moves along with the drop 
surface, see figure 2, we start from the continuity equation in axisymmetric (r,z) coordi¬ 
nates: 


= 0 . 


(Al) 


Ur dUr dUz 

r dr dz 

We now want to write the last two terms of the LHS of equation A1 in terms of the 
(||,_L)-coordinate system, that is: 


Ur dUr dUz Ur f dUr d _L dUr 5 || \ / dUz d _L dUz d II \ / A 

The two coordinate systems are related as follows (see also figure 2): 


11= —r sind -I- ZCOS0 
±= rcosO -I- zsind 

Using the relation above we can write equation A 2 as: 


(A3) 
(A 4) 


Ur ^ dUr ^ dUz Ur 
r dr dz r 


dUr „ dUr 


smi 


f duz 


+ 


\dH 


sind' 


duz 


COS0 


We now have to express Ur and Uz as function of (||,T), that is: 


(A 5) 


Ur{±, II) = ■u_l(-L, id cos6» - M||(T, id Sind (A6) 

W 2 (-L,II) = ^^(-L, |Dsind-h-«||(T, IDcosd (A7) 

Substituting the above expressions for Ur and Uz into equation A 5 and simplifying we 
find: 


Ur dUr dUz 
r dr dz 


Ur 

r 


du± ^ i9u|| 


(AS) 
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